#ifndef __MY_ANALYSIS_H__
#define __MY_ANALYSIS_H__

#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstdlib>
#include <map>
#include <algorithm>
#include <fstream>
#include <time.h>

#include "fastjet/ClusterSequence.hh"
#include "fastjet/tools/Pruner.hh"

#include "Njettiness.hh"

#include "CmdLine.hh"
#include "TFile.h"
#include "TChain.h"
#include "TProfile.h"

#include "TH1F.h"
#include "TH2F.h"

#include "Analysis.h"
#include "dipolarity.h"

using namespace std;
using namespace fastjet;

//used to sort jets via jet mass
class PseudoJet_mass_sorter
{
 public:
  bool operator()(const PseudoJet& j1, const PseudoJet& j2)
  {
    return j1.m() > j2.m();
  }
};


class MyAnalysis : public Analysis
{
 public:
  double evt_Mt;
  double evt_Ht;
  double evt_glu_m;
  TProfile* mass_profile;
  vector<Plot> jetmass_plots;
  
  MyAnalysis(const CmdLine& cmdline)
    : Analysis(cmdline)
    {
      evt_Mt=0;
      evt_Ht=0;
    }
 
  ~MyAnalysis(void)
    {
      cout<<"MESSAGE: "<<Event_+1<<" entries processed"<<endl;
      reset();
    }
 
  void hist_init();
  void hist_init2();
  void analyze();
  void fill_display();
  void fill_masspt();
  void fill_Njettiness();
  void fill_rutgers();
  void fill_rmass();
  void output();

  void analyze_glu(int);
  //determine whether event
  //passes our selections
  bool cut();
};


bool MyAnalysis::cut()
{
  //number of events before all cuts
  hist["nevt"]->Fill(0.0,weight);
  hist["nevt_unweighted"]->Fill(0.0);
  
  if(cmdline.present("-Mt"))
    {
      float Mt_cut=cmdline.value<float>("-Mt", 400);
      float Mt_present=0;
      for(int i=0; i<jets.size(); i++)
	if(jets[i].pt()>50)
	{Mt_present+=jets[i].m();}
      if(Mt_present < Mt_cut)
	return false;
    }

  if(cmdline.present("-jet4"))
    {
      if(Event_==0)
	{
	  cout<<"OPTIONS: Demand 1j > 100, at least 4j > 50"<<endl;
	  hist["nevt"]->GetXaxis()->SetBinLabel(2, "jet4");
	}

      if(jets.size()<4)
	return false;
      if(jets[0].pt()<100)
	return false;
      if(jets[3].pt() <50.0)
	return false;
    }

  if(cmdline.present("-jet2"))
    {
      if(Event_==0)
	{
	  cout<<"OPTIONS: Demand 1j > 100, at most 2j > 50"<<endl;
	  hist["nevt"]->GetXaxis()->SetBinLabel(2, "jet2");
	}

      if(jets.size()<2)
	return false;
      if(jets[0].pt()<100)
	return false;
      if(jets[1].pt()<50)
	return false;
      if(jets.size()>2 && jets[2].pt() >=50.0)
	return false;
    }

  if(cmdline.present("-jet2-tight"))
    {
      if(Event_==0)
	{
	  cout<<"OPTIONS: Demand 1j > 100, 1j > 50, at most 2j > 20"<<endl;
	  hist["nevt"]->GetXaxis()->SetBinLabel(2, "jet2-tight");
	}

      if(jets.size()<2)
	return false;
      if(jets[0].pt()<100)
	return false;
      if(jets[1].pt()<50)
	return false;
      if(jets.size()>2 && jets[2].pt() >=20.0)
	return false;
    }

  if(cmdline.present("-jet2-super"))
    {
      if(Event_==0)
	{
	  cout<<"OPTIONS: Demand 2j > 300, all other j < 20"<<endl;
	  hist["nevt"]->GetXaxis()->SetBinLabel(2, "jet2-supertight");
	}

      if(jets.size()<2)
	return false;
      if(jets[1].pt()<300)
	return false;
      if(jets.size()>2 && jets[2].pt() >=20.0)
	return false;
    }

  if(cmdline.present("-jet2-hard"))
    {
      if(Event_==0)
	{
	  cout<<"OPTIONS: Demand at least 2j > 300"<<endl;
	  hist["nevt"]->GetXaxis()->SetBinLabel(2, "jet2-hard");
	}

      if(jets.size()<2)
	return false;
      if(jets[1].pt()<300)
	return false;
    }

  
  hist["nevt"]->Fill(1.0,weight);
  hist["nevt_unweighted"]->Fill(1.0);
  return true;
}


void MyAnalysis::fill_display()
{
  if(i_key >= n_display)
    return;
  
  for(int i=0; i<InputJets.size(); i++)
    {
      //fill the eta phi energy deposit
      hist2[key_E[i_key]]->Fill(InputJets[i].eta(), 
				InputJets[i].phi());
      
      //fill the eta phi Pt desposit
      hist2[key_Pt[i_key]]->Fill(InputJets[i].eta(), 
				 InputJets[i].phi(),
				 InputJets[i].pt());
    }
  
  //loop over the first 4 jets
  //grab their constituents and make plots
  for(int i=0; i<4 && i<jets.size(); i++)
    {
      vector<PseudoJet> parts=jets[i].constituents();
      for(int k=0; k<parts.size(); ++k)
	hist2[key_jet[i_key][i]]->Fill(parts[k].eta(),
				       parts[k].phi(),
				       parts[k].pt());
    }
  
  
  analyze_glu(i_key);
  i_key++;
}

void MyAnalysis::fill_masspt()
{
  if(jets.size()>0){
    hist["jetpt1"]->Fill(jets[0].pt(),weight);
    hist["jetmass1"]->Fill(jets[0].m(),weight);
  }
  
  if(jets.size()>1){
    hist["jetpt2"]->Fill(jets[1].pt(),weight);
    hist["jetmass2"]->Fill(jets[1].m(),weight);
  }

  if(jets.size()>2){
    hist["jetpt3"]->Fill(jets[2].pt(),weight);
    hist["jetmass3"]->Fill(jets[2].m(),weight);
  }
  
  double Mt=0, Ht=0;
  int njet=0;
  
  for(int i=0; i<jets.size(); i++)
    {
      hist["jetpt"]->Fill(jets[i].pt(),weight);
      hist["jetmass"]->Fill(jets[i].m(),weight);
      if(jets[i].pt()>50)
	{
	  njet++;
	  Mt+=jets[i].m();
	  Ht+=jets[i].mperp();
	}
    }
  hist["njet"]->Fill(njet,weight);
  hist["Mt"]->Fill(Mt,weight);
  hist["Ht"]->Fill(Ht,weight);
  hist2["Mt_Ht"]->Fill(Ht, Mt,weight);

  evt_Mt=Mt;
  evt_Ht=Ht;
}

void MyAnalysis::fill_rmass()
{
  for(int i=0; i<jetmass_plots.size(); i++)
    {
      makejets(antikt_algorithm, i*0.2+0.6);
      double M_J = 0;
      for(int j=0; j<jets.size(); j++)
	{
	  if(jets[j].pt()<50.0)
	    break;
	  M_J+=jets[j].m();
	}
      jetmass_plots[i]->Fill(M_J, weight);
      mass_profile->Fill(i*0.2+0.6+0.001, M_J, weight);
      hist2["Mt_radius"]->Fill(i*0.2+0.6+0.001, M_J, weight);
    }
}



void MyAnalysis::fill_Njettiness()
{
  //setup N-subjettiness   
  //power for the angular dependence
  //2 = quadratic
  double beta=1.0;
  double R0=1.2; 
  //R0 used only for normalization;
  double Rcut=100;
  
  NsubParameters param(beta, R0, Rcut);
  Njettiness nSubMin(Njettiness::min_axes,
		     param);


      //fill Nsubjettiness for first 2 jets
      if(jets.size()>0)
	{  
	  //get tau for hardest jet
	  vector<PseudoJet> jet1_sub=jets[0].constituents();
	  double jet1_tau1=nSubMin.getTau(1,jet1_sub);
	  double jet1_tau2=nSubMin.getTau(2,jet1_sub);
	  double jet1_tau3=nSubMin.getTau(3,jet1_sub);
	  
	  hist["jet1_tau1"]->Fill(jet1_tau1,weight);
	  hist["jet1_tau2"]->Fill(jet1_tau2,weight);
	  hist["jet1_tau3"]->Fill(jet1_tau3,weight);
	  
	  double c_tau21=jet1_tau2/jet1_tau1;
	  double c_tau32=jet1_tau3/jet1_tau2;
	  double c_tau31=jet1_tau3/jet1_tau1;

	  double c_pt=jets[0].pt();
	  double c_m=jets[0].m();
  
	  hist["jet1_tau21"]->Fill(c_tau21,weight);
	  hist["jet1_tau32"]->Fill(c_tau32,weight);
	  hist["jet1_tau31"]->Fill(c_tau31,weight);
	  
	  hist2["jet1_tau3_tau2"]->Fill(jet1_tau2, jet1_tau3,weight);
	  
	  hist2["jet1_mass_tau21"]->Fill(c_tau21,
					 c_m, weight);
	  hist2["jet1_mass_tau31"]->Fill(c_tau31,
					 c_m, weight);
	  hist2["jet1_mass_tau32"]->Fill(c_tau32,
					 c_m, weight);
	  
	  hist2["jet1_pt_tau21"]->Fill(c_tau21,
				       c_pt, weight);
	  
	  hist2["jet1_pt_tau31"]->Fill(c_tau31,
				       c_pt, weight);
	  
	  hist2["jet1_pt_tau32"]->Fill(c_tau32,
				       c_pt, weight);
	  

	  
	  

	  if(c_pt <200.0)
	    {
	      hist["jet1_tau31_ptbin1"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_ptbin1"]->Fill(c_tau21, weight);
	    }
	  else if(c_pt <600.0)
	    {
	      hist["jet1_tau31_ptbin2"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_ptbin2"]->Fill(c_tau21, weight);
	    }
	  else if(c_pt <1200.0)
	    {
	      hist["jet1_tau31_ptbin3"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_ptbin3"]->Fill(c_tau21, weight);
	    }
	  else if(c_pt <1800.0)
	    {
	      hist["jet1_tau31_ptbin4"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_ptbin4"]->Fill(c_tau21, weight);
	    }
	  else
	    {
	      hist["jet1_tau31_ptbin5"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_ptbin5"]->Fill(c_tau21, weight);
	    }


	  if(c_m > 50 && c_m <150)
	    {
	      hist["jet1_tau31_mbin1"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_mbin1"]->Fill(c_tau21, weight);
	    }
	  
	  if(c_m > 145 && c_m <205)
	    {
	      hist["jet1_tau31_mbin2"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_mbin2"]->Fill(c_tau21, weight);
	    }
	  
	  if(c_m > 200 && c_m <400)
	    {
	      hist["jet1_tau31_mbin3"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_mbin3"]->Fill(c_tau21, weight);
	    }

	  if(c_m >400)
	    {
	      hist["jet1_tau31_mbin4"]->Fill(c_tau31, weight);
	      hist["jet1_tau21_mbin4"]->Fill(c_tau21, weight);
	    }

	  hist2["Mt_tau21"]->Fill(c_tau21, evt_Mt, weight);
	  hist2["Mt_tau32"]->Fill(c_tau32, evt_Mt, weight);
	  hist2["Mt_tau31"]->Fill(c_tau31, evt_Mt, weight);



	}

      if(jets.size()>1)
	{
	  vector<PseudoJet> jet2_sub=jets[1].constituents();
	  double jet2_tau1=nSubMin.getTau(1,jet2_sub);
	  double jet2_tau2=nSubMin.getTau(2,jet2_sub);
	  double jet2_tau3=nSubMin.getTau(3,jet2_sub);
	  
	  hist["jet2_tau1"]->Fill(jet2_tau1,weight);
	  hist["jet2_tau2"]->Fill(jet2_tau2,weight);
	  hist["jet2_tau3"]->Fill(jet2_tau3,weight);
	  
	  hist["jet2_tau21"]->Fill(jet2_tau2/jet2_tau1,weight);
	  hist["jet2_tau32"]->Fill(jet2_tau3/jet2_tau2,weight);
	  hist["jet2_tau31"]->Fill(jet2_tau3/jet2_tau1,weight);
	  
	  hist2["jet2_tau3_tau2"]->Fill(jet2_tau2, jet2_tau3,weight);
	}


      //grab the most massive jet
      if(jets.size()>0)
	{
	  vector<PseudoJet> jets_m=jets;
	  sort(jets_m.begin(), jets_m.end(), PseudoJet_mass_sorter());
	  vector<PseudoJet> jet1_sub=jets_m[0].constituents();
	  double jet1_tau1=nSubMin.getTau(1,jet1_sub);
	  double jet1_tau2=nSubMin.getTau(2,jet1_sub);
	  double jet1_tau3=nSubMin.getTau(3,jet1_sub);
	  
	  double c_tau21=jet1_tau2/jet1_tau1;
	  double c_tau32=jet1_tau3/jet1_tau2;
	  double c_tau31=jet1_tau3/jet1_tau1;


	  hist2["Mt_m_tau21"]->Fill(c_tau21, evt_Mt, weight);
	  hist2["Mt_m_tau31"]->Fill(c_tau31, evt_Mt, weight);
	    
	}
      
      
}





void MyAnalysis::fill_rutgers()
{
  makejets(antikt_algorithm, 0.5);
  double min_pt_cut=30.0;
  int index=0;
  for(; index<jets.size(); index++)
    if(jets[index].pt() < min_pt_cut)
      break;

  //6 jets cut
  if(index < 6) return;

  double mass=-1;
  double sumpt=-1;
  double diag=-1;

  for(int i=0; i<index; i++)
    for(int j=0; j<i; j++)
      for(int k=0; k<j; k++)
	{
	  double c_mass=(jets[i]+jets[j]+jets[k]).m();
	  double c_sumpt=jets[i].pt()+jets[j].pt()+jets[k].pt();
	  hist2["Mjjj_sumpt"]->Fill(c_sumpt, c_mass, weight);

	  if(c_sumpt - c_mass > 100)
	    {
	      hist["Mjjj_100"]->Fill(c_mass, weight);
	      if(c_sumpt - c_mass > 200)
		{
		  hist["Mjjj_200"]->Fill(c_mass, weight);
	      	  if(c_sumpt - c_mass > 300)
		    {
		      hist["Mjjj_300"]->Fill(c_mass, weight);
		      if(c_sumpt - c_mass > 400)
			{
			  hist["Mjjj_400"]->Fill(c_mass, weight);
			  if(c_sumpt - c_mass > 500)
			    hist["Mjjj_500"]->Fill(c_mass, weight);
			}
		    }
		}
	    }
	  
		      


	  if(c_sumpt - c_mass > diag)
	    {
	      diag= c_mass - c_sumpt;
	      mass=c_mass;
	      sumpt=c_sumpt;
	    }
	}
  hist2["Mjjj_sumpt_max"]->Fill(sumpt, mass, weight);
}


void MyAnalysis::output()
{
  

  double beta=1.0;
  double R0=1.2; 
  //R0 used only for normalization;
  double Rcut=100;
  
  NsubParameters param(beta, R0, Rcut);
  Njettiness nSubMin(Njettiness::min_axes,
		     param);


  for(int i=0; i<jets.size(); i++)
    {
      cout<<"jets "<<i<<": "
	  <<jets[i].px()<<", "
	  <<jets[i].py()<<", "
	  <<jets[i].pz()<<", "
	  <<jets[i].e()<<", ";
	  


	  //get tau for hardest jet
	  vector<PseudoJet> jet1_sub=jets[i].constituents();
	  double jet1_tau1=nSubMin.getTau(1,jet1_sub);
	  double jet1_tau2=nSubMin.getTau(2,jet1_sub);
	  double jet1_tau3=nSubMin.getTau(3,jet1_sub);

	  cout<<jet1_tau1<<", "
	      <<jet1_tau2<<", "
	      <<jet1_tau3<<", "<<endl;
	    
    }
  
}

#endif

